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ABSTRACT 



An inverse model involving AVHRR imagery and the heat equation with dynamical 
constraints on the divergence, kinetic energy and vorticity of the solutions was used by 
Kelly (1989) to produce velocity fields that were in good agreement with Acoustic 
Doppler Current Profiler (ADCP) data. Dynamic heights derived from GEOSAT radar 
altimeter data have also been used to determine near-surface geostrophic currents. 
Synthetic GEOSAT-derived velocity data w r as generated from ADCP data collected as 
part of the Coastal Transition Zone (CTZ) Field Program. The inverse model was run 
with AVHRR imagery that w r as coincident to the CTZ Field Program ADCP data and 
the synthetic velocity data was added as an additional constraint on the model's sol- 
ution. The resulting velocity solutions w r ere much improved over those given by the in- 
verse model alone. Refinement of this method involving a combination of different data 
sources should improve efforts to determine near-surface velocities of the ocean entirely 
by remote means. 
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I. BACKGROUND 



A. FEATURE TRACKING 

From the first availability of sequential satellite images, researchers have been mo- 
tivated towards estimating velocities of the atmosphere and ocean on a large scale from 
the motion of features in the images. The first attempts at estimating near-surface ve- 
locity fields in the ocean were by means of a technique called feature tracking. 

Early feature tracking involved manually tracking the movement of patterns that 
could be easily identified in consecutive images. The displacement of the features be- 
tween consecutive images was used to determine velocities. The identification and 
tracking of clouds in visible or infrared imagery has been used to calculate winds in the 
atmosphere. The movement of chlorophyll concentrations in Coastal Zone Color 
Scanner (CZCS) imagery and Sea Surface Temperature (SST) gradients in infrared im- 
agery have both been used to estimate near-surface current velocities in the ocean. 

Svejkosky (1988) found RMS error of 0.06 ms -1 between feature tracking derived 
velocities from Advanced Very High Resolution Radiometer (AVHRR) imagery and 
those obtained from sea surface drifters. The error was low because constraints on the 
feature tracking were restrictive, allowing use of only the most distinctive features and 
holding the temporal variation between the imagery and surface drifters to five hours. 

There are several drawbacks to manual feature tracking. It is very time consuming 
and highly subjective. Different operators will achieve differing results from the same 
sets of imagery. Feature tracking cannot be used on low contrast or low gradient images 
since they have no features discernable to the human eye. 

Feature tracking has been made more objective through computer automation. A 
technique developed by Leese et al. (1971) for cloud tracking was adapted for use with 
SST images by Emery et al. (1986). The Maximum Cross Correlation (MCC) technique 
involves cross correlating small sections in the first image with slightly larger sections in 
the second image. The cross correlated "windows" in each image have the same ge- 
ographical center point. Displacement vectors from the window centers to their respec- 
tive points of maximum correlation represent the surface velocity field. 

The main drawback of the MCC technique is that it is adversely affected by the 
distortion and rotation of the features as they move. Additionally, the technique is un- 
able to distinguish motion that is occurring parallel to the thermal or color gradients. 
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Tokmakian et al. (1990) achieved RMS errors of 0.18-0.23 ms -1 using the MCC 
technique with AVHRR images having 12 hour separations. They were able to better 
quantify the actual error in the MCC technique by applying the MCC technique to 
synthetic SST images. Synthetic SST images were produced by advecting the temper- 
ature field of an infrared satellite image with the velocity field of a numerical model and 
then sampling the model output at the desired intervals. The error of the MCC tech- 
nique could then be found directly by comparing its velocities to the known velocities 
of the numerical model. Analysis indicated that the error might be reduced by cutting 
the image separation to 6 hours and by incorporating some means to account for rota- 
tion. They also felt that improvement might be realized by starting the correlative 
searches in high gradient regions and by somehow eliminating random, spurious high 
correlations that were occurring between distant, unrelated areas of the images. 

B. INVERSE METHODS 

A different objective approach to deriving near-surface velocities from AVHRR im- 
agery is through the use of inverse methods. Wunsch (1978, 1985), Fiadeiro and 
Veronis (1982, 1984) and Roemmich (1979, 1981) are a few of those who have applied 
inverse methods to measured distributions of tracers in the ocean in order to infer large 
scale circulation and current velocities. 

Generally speaking, inverse problems involve combining measurements of conserved 
properties and their corresponding governing equations into a mathematical system. 
The resulting system can be underdetermined or overdetermined and is solved using nu- 
merical methods for a velocity field which could have produced the measured distrib- 
ution of properties. 

A critical assumption in the subjective and objective feature tracking techniques 
described earlier was that the temperature gradients were being advected without signif- 
icant change. This is actually not the case as surface heat fluxes, heat diffusion and 
vertical mixing of the water column all act to change the surface temperatures and gra- 
dients. Inverse methods would seem to be better suited to this problem because they 
better represent the actual physical processes involved. Source terms as well as other 
physics can be included in the model rather than assuming straight advection. 

The relevant forward problem involving temperature advection and the heat 
equation is to calculate some future surface temperature field given an initial surface 
temperature field and a horizontally advecting velocity field. The corresponding inverse 
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problem is to determine the advecting velocity field given only the initial and final tem- 
perature fields. Even the forward problem in this case is not well understood since in 
reality it is a three-dimensional non-conservative problem. Kelly (1989) applied an in- 
verse model to the heat equation and found good agreement between the model's ve- 
locity solutions and coincident Doppler Acoustic Log (DAL) data.l More detail 
concerning the inverse method with temperature advection and the heat equation will 
be presented in Chapter II. 

C. RADAR ALTIMETER DERIVED NEAR-SURFACE VELOCITIES 

An entirely different means of deriving near-surface velocities from satellite data is 
to use radar altimeter data. Radar altimeters accurately measure the distance from the 
satellite to the surface of the earth. Sharp pulses of electromagnetic radiation (at a fre- 
quency of 13.5 Ghz for the GEOSAT radar altimeter) are beamed towards the earth's 
surface and their travel time converted to distance. Corrections described in Cheney et 
al. (1987) for the GEOSAT altimeter are applied to the raw data to correct for tides, 
atmospheric water vapor, tropospheric and ionospheric delays and surface atmospheric 
pressure. Over the ocean, subtracting the satellite radar altimeter measurement from the 
height of the satellite's orbit results in the height of the sea surface relative to the center 
of the earth. Figure 1 illustrates these relationships. The sea surface height variability 
is due mainly to a combination of 1) the shape of the geoid, 2) ocean tides, and 
3) ocean currents. Having already corrected for tidal effects, a correction must be ap- 
plied for the variations due to the geoid in order to realize the sea surface height vari- 
ations resulting solely from the ocean's currents. 

The true geoid is the gravitational equipotential surface corresponding to the level 
the ocean would seek with gravity and centrifugal forces as the only forces acting upon 
it. Geoid variations range from —104 m to +60 m while the variations in dynamic 
height due to ocean currents are generally +1.5 m. Averaging the radar altimeter 
measurements of many orbits over the same location will determine an altimetric geoid. 
An altimetric geoid is often used as an approximation for the true geoid. Unfortunately, 
such time averaging includes sea surface height variations due to the mean flow'' of the 
ocean in the altimetric geoid. Other means can be used to find the true geoid over the 
ocean. The long wave components of the geoid can be found through their effect on 



1 The DAL is an earlier version of the Acoustic Doppler Current Profiler (ADCP) which will 
be described later. Both will be referred to as ADCP in this paper. 
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satellite orbits. The short wave components can be found through detailed gravity sur- 
veys. Once again difficulty arises because sufficiently detailed gravity surveys are avail- 
able for only small areas of the world's oceans. Additionally, security classification 
restrictions have precluded the use of an accurate geoid in general scientific research 
(Stewart, 1985). 

Roemmich and Wunsch (1982) describe +10 cm as an attainable accuracy for a 
geoid obtained by using inverse methods with radar altimeter and hydrographic data to 
remove the mean oceanographic signal from an altimetric geoid. They also demonstrate 
that 10 cm is the maximum allowable error for a true geoid in order for it to be useful 
in improving estimates of the general circulation of the ocean. 

Joyce et al. (1986) project a capability for significant improvement in the geoid 
through the use of inverse methods with a combination of hydrographic and ADCP ve- 
locity data. Near-surface geostrophic velocities obtained from the inversion would be 
used to calculate the dynamic height of the ocean's surface. The resulting dynamic 
heights would be used to further refine the geoid to accuracies not attainable using 
gravimetry. 

The dynamic height variation of the ocean that remains after the final geoid cor- 
rection is applied to the radar altimeter data. The dynamic height variation along the 
satellite track bears the following relationship to the cross-track geostrophic current ve- 
locity 



~g dh 

f dy 



(l.i) 



where u x is the cross-track surface velocity, g is gravity, /is the Coriolis parameter and 
is the dynamic height variation along the satellite track. The accuracy of radar 
altimeter derived geostrophic velocities is related to the accuracy of the geoid through 
its effect upon the accuracy of the dynamic height field. 

Kelly and Gille (1990) used information from the GEOSAT radar altimeter to cal- 
culate Gulf Stream velocities that were in good agreement with shipboard ADCP meas- 
urements. They subtracted the mean of all of the height profiles as a substitute for a 
reference geoid in order to get GEOSAT residual height profiles. A theoretical model 
of the Gulf Stream velocity profile was used to generate synthetic height profiles. The 
synthetic height profiles were then added to the GEOSAT residual height profiles in or- 
der to account for the mean flow of the Gulf Stream. They achieved an RMS error of 
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0.23 ms -1 between a GEOSAT derived velocity profile and one that had been measured 
with an ADCP four days earlier. Most of the error was attributed to the temporal var- 
iability of the Gulf Stream evidenced through sequential ADCP profiles. 

D. RADAR ALTIMETER DERIVED VELOCITIES AS CONSTRANTS ON AN 
INVERSE MODEL 

The system of equations and unknowns comprising an inverse model may be sup- 
plemented by known data. Such a system is said to be constrained by the known data 
which acts to guide the inverse model towards a solution that is representative of the 
known data. This project will evaluate the effects of adding simulated radar altimeter 
derived velocity data as an additional constraint to the inversion. 
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II. METHODS 



Kelly (1989) applied an inverse model using AVHRR imagery as input to the heat 
equation 



r f + u.vr=K-v 2 r+5 



( 2 . 1 ) 



where T is temperature, u is velocity, k is a diffusion coefficient, 5 is a source term and 
subscripts denote partial differentiation. 

The horizontal temperature gradients were computed from 256x256 pixel images by 
a linear least-squares fit across 16x16 pixel tiles. With an AVHRR sensor footprint of 
1.15 km the resulting a 18x18 km boxes could sustain a maximum cross-isotherm ve- 
locity of 15 kmd-‘ over a one day time separation without violating the Courant- 
Friedrich-Lewev (CFL) criterion Ax > uAt. The term T, was approximated by the finite 

difference — for sequential images separated in time. A small correction for bias 

in the images was made by calculating a mean image temperature difference using only 
low gradient tiles. That mean was then subtracted from the temperature difference cal- 
culated for each tile. In other words, the low gradient regions were assumed to be con- 
stant temperature and any temperature change occurring in those regions over the image 
interval was assumed to be bias in the images. The source term 5 accounted for surface 
heat fluxes and vertical entrainment which were assumed to vary slowly over the image. 
It also included large scale temperature error introduced into the satellite imagery by 
atmospheric effects. Diffusion was neglected as being at least an order of magnitude 
smaller than advection. 

Considering u as across and along isotherm components showed that the advective 
temperature change was due to the cross isotherm component only. 



u .vr=(u f + u a ).vr=u,.vr 



( 2 . 2 ) 



Since there was no contribution to advection by u„ (or either component where there 
were negligible temperature gradients) an infinite number of solutions to 2.1 existed. In 
order to achieve solutions that were unique but also physically reasonable, limitations 
on the divergence, energy and vorticity of the inversion solution were chosen as 
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constraints on the problem. The equation to be minimized was (Kelly, 1989) 

JX I 7) + u • V7~ — S ly + a 2 | V • u 1^/ + /? 2 | u + y 2 | V x u |^] (2.3) 

<7 

where T, + u^VT — S is the heat equation with diffusion neglected, V • u is the 
divergence, u 2 is the energy and V x u is the vorticity of the solution. Weighting factors 
a, p and y are used to specify the relative importance of each constraint to the inversion. 

Minimization of divergence forced temperature changes to be made through move- 
ment of features rather than by convergence/divergence. Minimizing the energy of the 
solution most closely approached the solution allowing only cross-isotherm flow. A 
minimum vorticity solution allowed jets only where necessary to satisfy the heat 
equation. 

The resulting overdetermined system of the heat equation and constraints was 
solved using a matrix factorization algorithm. The algorithm used least squares meth- 
odology in order to satisfy all equations as equally as possible. 

Testing of the inverse model was done using synthetic images with synthetic velocity 
fields and then real images with synthetic velocity fields. The purpose of the tests was 
to find the proper weightings for the divergence and energy constraints as well as deter- 
mine the model's sensitivity to some types of expected error. 

For the first test, a relatively low gradient (< 0. l°Kknr‘) synthetic temperature im- 
age was advected by a variety of uniform or slightly divergent synthetic velocity fields 
( V max 15 kmd~ l ). The temperature differences of the resulting images were then 
input to the model. Varying the weights on the constraints revealed their effects on the 
inverse model in terms of percent misfit of the model's solution to the actual time-rate 
of temperature change. Analysis indicated that the main source of error when utilizing 
"perfect" data was underparameterization of the velocity field by the four term Fourier 
expansion. Solution velocities were also lower than the actual velocities due to the low 
gradients of the synthetic temperature image. Negligible temperature changes between 
the images did not require advecting velocities in order to solve the heat equation. 

The second test involved the advection of a real temperature image with a divergent 
synthetic velocity field. The best solutions from that test had lower divergence and 
energy than the actual solution. There was also a tendency for the model to develop 
spurious structure in the solutions. The velocity error was smaller however for the 
solutions in the second test sequence because the real image had stronger gradients than 
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the synthetic image used in the first set of tests. Noise calculated on the basis of imagery- 
calibration error analysis was then added to the test data. Solutions again hin!d lower 
divergence and energy than the actual solution. 

For inversions with real data, weights on the constraints were chosen to yield a 
family of solutions within a certain range of temperature misfit values. The acceptable 
range of misfit values was determined from estimates of the image calibration error and 
the small-scale temperature error estimated from the model's large-scale source terms 
S. Inversions performed on real images produced solutions that were for the most part 
qualitatively correct with respect to coincident ADCP data but having lower velocities, 
divergence and energy. 

Additional weighted constraints in the form of known velocities or estimated veloc- 
ities from another data source could also be added to the system. If the known or esti- 
mated velocities happened to be representative of the actual along isotherm velocities, 
significant improvement should be realized over the inversion alone. 
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III. DATA 



A. AVHRR IMAGERY 

The AVHRR images used for this project were four clear, consecutive images from 
the NOAA 11 satellite. They were collected 12 hours apart over the period 
16-17 July 1988 (Julian dates 198 and 199). Images from the NOAA 10 satellite were 
also available for this period but were not used because the NOAA 10 AVHRR lacked 
the Channel 5 data necessary to use the Multi-channel Sea Surface Temperature 
(MCSST) algorithms. Table 1 provides information pertinent to the imagery. The im- 
ages can be seen as the backgrounds for the Coastal Transition Zone (CTZ) ADCP ve- 
locity field (Figure 6) and the inversion-alone solutions (Figures 7-9). 



Table I. NOAA 11 AVHRR IMAGERY 







Time of Pass 


Satellite 

Zenith 

Angle 


CTZ Temperature Statistics 


Image 

Designation 


Julian 

Date 


GMT 


Local 


T maX 


T min 


T avg 










(°C) 


(°C) 


(°C) 


19812 


198 


1200 


0500 


35° 


17.8 


7.1 


14.6 


19823 


19S 


2300 


1600 


3° 


17.4 


7.5 


13.9 


19912 


199 


1200 


0500 


45° 


18.9 


7.0 


13.9 


19923 


199 


2300 


1600 


18° 


18.4 


7.9 


13.7 



The images were navigated by Mark Abbott at Oregon State University to an ac- 
curacy of ±1 pixel using Global Imaging, Inc. satellite data processing software. The 
current National Environmental Satellite, Data and Information Services (NESDIS) 
operational MCSST algorithms for the NOAA 11 satellite were used to calculate sea 
surface temperatures from the AVHRR Channel 4 and 5 radiances. Values outside the 
range 7-19°C were masked and rejected as not being representative sea surface temper- 
atures. That range was based on data from Huyer et al. (1990) for the Coastal Transi- 
tion Zone (CTZ) during this time, expanded to account for the AVHRR RMS 
temperature error of 0.5°C reported by McClain et al. (1985). 

No cloud masking was done on the images because of their high clarity. There was 
an area of clouds in Image 19812 that covered a small area of the extreme southwest 
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comer of the model grid. The clouds did result in erroneous velocities from the inverse 
model in that region. There also appeared to be coastal fog present in Images 19823, 
19912 and 19923. The presence of fog would account for the abnormally high nearshore 
velocities in some of the model solutions using those images. 

B. CTZ DATA 

The control data for this project consisted of Acoustic Doppler Current Profiler 
(ADCP) and Conductivity, Temperature and Depth (CTD) data collected for the 
Coastal Transition Zone (CTZ) Field program over the period 13-18 July 1988 by the 
R/V PT SUR. Table 2 summarizes the temporal aspects of the CTZ88 Grid II ADCP 
data, collected in seven northwest-southeast transects labelled A-G. 



Table 2. CTZ88 GRID II ADCP SURVEY 



Line 


Date 

(1988) 


Time (GMT) 


Start 


End 


A 


13-14 July 


0813 


0201 


B 


14 July 


0446 


1945 


C 


14-15 July 


2211 


1015 


D 


15-16 July 


1239 


0215 


E 


16-17 July 


0504 


0124 


F 


17 July 


0417 


2000 


G 


17-18 July 


2346 


1515 



The ADCP velocity field and dynamic height fields of the CTZ during this period (Fig- 
ures 2 and 3) show good agreement. The dominant feature in both data sets is a strong 
surface jet extending northeast to southwest across the study region, with maximum 
speeds reaching =* 90 cms -1 . The circulation in the northwest comer is weakly cyclonic. 
A weaker onshore "return" flow exists to the south of the main jet. 

An ADCP measures current velocity at various depths (bins) by range-gating a so- 
nar pulse and recording the Doppler shift of the signal returned from each bin 
(Figure 4). Calculation and subtraction of the effects of the ship's motion from the total 
Doppler shift velocity leaves the velocity due to the motion of scatterers in the ocean. 
The scatterers are presumed to be drifting with the current. Averaging over some opti- 
mal length of time yields a mean current for each depth bin over an area corresponding 
to the distance travelled by the ship over the averaging time. Error in the ADCP velocity 
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is mainly due to signal processing difficulties arising from the orders of magnitude dif- 
ference that can exist between the actual currents and the ship's speed. Other error re- 
sults from the ship's motions due to wave action and the detrimental effects of 
navigation error on calculations of the ship's speed (Kosro, 1985). 

C. SYNTHETIC GEOSAT-DERIVED VELOCITY DATA 

Synthetic GEOSAT-derived velocity data, generated from ADCP control data, was 
used as an additional constraint on the inversion. The cross-track ADCP data were 
averaged over a time interval that approximated the 7.5 km sensor footprint of 
GEOSAT. The ADCP survey lines were used to simulate GEOSAT ground tracks. The 
line spacing of the ADCP survey was 25 km compared to the 2 s 130 km track spacing 
of repeat GEOSAT orbits over the CTZ. The angle between the ADCP survey lines and 
the model grid was 30° while the angle between the GEOSAT ground tracks and the 
model grid would have been s£ 25°. Figure 5 shows the synthetic GEOSAT-derived ve- 
locity field. 

The error of the ADCP data set was estimated at 5-10 cms -1 by Huyer et al. (1990). 
Since the synthetic velocity data was derived from actual velocity data rather than syn- 
thetic raw altimeter data or synthetic dynamic height profiles, any error normally in- 
curred in the process of converting raw altimeter data to geostrophic velocities was 
avoided. The availability of an accurate geoid was not a factor and no estimate of the 
geoid was made. Although the ADCP data was collected over time and therefore was 
not coincident to all of the image pairs used with the model, for the purposes of this 
study, the synthetic velocity data were assumed to be coincident with the imagery'. 
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IV. THE MODEL 



The model used was essentially the same one used and described in Kelly (1989). 
The model had been previously modified to accept known velocities in the form of 
ADCP data which was weighted and used as a constraint on the inversion. The newer 
model lacked the vorticity constraint of the earlier model which was removed because it 
had only a minor effect on the results. Additionally, the newer model parameterized the 
velocities as biharmonic splines rather than as Fourier coefficients. The use of 
biharmonic splines reduced the tendency of the solution to resonate around areas where 
data were masked (Kelly, personal communication). 

In order to incorporate the synthetic GEOSAT-derived velocity data as a constraint 
on the inversion, the model was further modified to minimize the difference between the 
velocity component normal to the satellite track and the inversion velocity. The inver- 
sion velocities were free to vary in the direction of the unknown component along the 
satellite track. The family of equations described in 2.1 became 

Xf I T, + u • VT - S \jj + a 2 1 V . u fij + P 2 1 u | ]j + t ] 2 1 v, - (u cos 6 + v sin 6) (4. 1) 

‘J 

where the terms are as described in 2.1 with the addition of another weighting function 
(rj), the synthetic GEOSAT cross-track velocity component (v,) and the angle between 
the model's grid and the satellite track (6). 
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V. RESULTS 



A. THE INVERSE MODEL WITHOUT VELOCITY CONSTRAINTS 

The modified model was run without velocity constraints. The resulting velocity 
fields were products of the inversion alone with constraints only on the divergence and 
energy of the solution. The model was run on each of the three image pairs with 
12 hour separations and also on the two image pairs with 24 hour time differences. 

The weights on the divergence and energy of the solution (a 2 and / ? 2 in equation 2.3) 
were held constant at .01 and .03 for all model runs. These values were used because 
they gave qualitatively good results and project time constraints precluded a search for 
weights that might have yielded better solutions. Better results (quantitative reductions 
in the velocity error between the model output and known solutions) may be possible 
through further optimization of a and /?. 

Several characteristics were calculated for each of the velocity solutions in order to 
allow quantitative comparisons. The energy and divergence of each velocity solution 
were computed. A percentage temperature error for each velocity solution was com- 
puted using the following equation: 



1 err 




I 7~, — ((u • VT) + S) \jj 
I T, \jj 



(5.1) 



where 7) is the temperature difference between the image pairs, u is the model's velocity 
solution, VT is the computed image gradient, S is the source terms calculated by the 
model and i=j= 16. A percentage velocity error between the model velocity solution 
and the CTZ88 Grid II ADCP velocity field was calculated in the following manner 




(“« ~ “J 2 + ( v « ~ v m) 2 

K + v mf 



(5.2) 



where n ( = 75) was the number of tiles where measured velocity data existed and sub- 
scripts a and m denote ADCP and model velocities respectively. Table 3 provides a 
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summary of the characteristics of the solutions using the inversion without any known 
velocities as constraints. Figure 6 is the CTZ88 Grid II ADCP velocity field overlaid 
on Image 19823. Figures 7-11 are the velocity fields from the inversion alone overlaid 
on either the second image of the pair used, or the image at the midpoint of the interval 
between the image pair. 



Table 3. CHARACTERISTICS OF INVERSION-ALONE SOLUTIONS 



Image pair 


Energy 

(km 2 d 2 ) 


Divergence 

(d- 1 ) 


T.„ (%) 


v„ (%) 


19812-19823 


180.0 


0.316 


5.94 


92.4 


19823-19912 


259.0 


0.421 


34.20 


107.0 


19912-19923 


293.0 


0.579 


7.76 


101.0 


19812-19912 


107.0 


0.314 


11.20 


97.9 


19823-19923 


96.6 


0.276 


7.39 


96.0 



The strong jet from the CTZ88 Grid II ADCP and dynamic height fields 
(Figures 1 and 2) is also evident in the AVHRR imagery. The cyclonic eddy is also 
apparent to the west of the cold front near Point Arena. The eddy is partially included 
in the northeast quadrant of the model grid although it is not well sampled by the 
CTZ88 Grid II ADCP survey. In all cases, the inverse model when used alone appeared 
to have problems handling the jet. The solutions do show some of the convergence 
necessary' for the formation of the jet, but don't actually define it. Part of the difficulty 
may have been due to the narrowness of the cold filament in relation to the scale of the 
model grid and because the jet velocities were predominately along-isotherm. 

Looping the images revealed movement that was not readily discemable from the 
still pictures. Over the sequence of images, the middle north-south oriented portion of 
the cold front translated westward =* 18 km, a distance roughly equal to its width. The 
east-west oriented region of the cold filament remained nearly stationary' and appeared 
to narrow and get colder over the image sequence. The looped sequence also showed 
the growth and southeastward movement of the northern cylconic eddy. The movement 
of the filament over the interval between images appeared to have weakened the average 
gradients computed by the model causing it to inadequately define the jet in that area, 
producing the velocity "holes” in the center of most solutions. In other words, the 
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inversion-alone vectors in many cases were simply an accurate representation of the 
cross isotherm velocities (translation of features). 

The inversion-alone solution for image pair 19812-19823 (Figure 7) gave the best 
quantitative results in terms of lowest temperature and velocity error (Table 3). It also 
looked the best qualitatively because of its lack of erroneous structure in the low' gradi- 
ent regions. The solution did have the "hole" associated with the narrow east-west sec- 
tion of the jet. Erroneous velocities were present in the extreme southwestern corner of 
the model grid owing to the cloud contamination in Image 19812. The possible presence 
of coastal fog could also call into question some of the nearshore velocities. 

The solution for image pair 19823-19912 (Figure 8) had the largest temperature and 
velocity errors of the inversion-alone solutions. This image pair had the largest change 
in T m „ between any of the pairs as well as the second largest change in T mjn (Table 1). 
This temperature change is the opposite of what would be expected from any diurnal 
effect resulting from the 12 hour day/night image separation. This image pair also has 
the largest satellite zenith angle difference (Table 1) of any of the pairs. The zenith angle 
difference or alternatively just the large satellite zenith angle of Image 19912 could be 
causing a problem through the MCSST algorithm used in these cases. The a* 9 km 
westward translation of the north-south section of the jet between the images was the 
largest of any of the pairs. The area of the bend from north-south flow to east-west flow 
also translated southeastward. These feature translations could be responsible for the 
strong eastward (cross-jet) velocities over the north-south segment of the jet and the 
seemingly erroneous southward velocities in the low gradient region to the north and 
west of the bend in the jet. The noise present at the bottom of Image 19912 was rejected 
as out-of-range temperatures and did not have any effect on the model. 

Image pair 19912-19923 produced a solution (Figure 9) that fell between the pre- 
ceding two descriptions. Feature translation did not appear to occur between the images 
in this instance but there were rather severe problems with erroneous velocities in the 
relatively low gradient region to the west of the jet. 

Feature translation appeared to be the largest error source in the solutions using 
24 hour separation (Figures 10 and 11). In these cases the greater image temporal 
separation allowed the features to translate further over the interval between images. 
As expected considering the previous results, the solution involving Image 19912 
(Figure 1 1) had the worst characteristics. Although no diurnal effects were apparent in 
any of these cases, the 24 hour image separation would minimize any adverse effects 
that a diurnal temperature variation would have on the inverse model. 



15 



B. THE INVERSE MODEL WITH SYNTHETIC VELOCITY CONSTRAINTS 

The next step was to add the synthetic GEOSAT-derived velocity data as an addi- 
tional constraint on the inversion. The modified model was run and various combina- 
tions of the synthetic velocity data lines B-G were added to the problem. The actual 
energy of the CTZS8 Grid II ADCP velocity field was calculated to be ^ 875 km 2 d~ J , 
however velocity solutions from the model with this energy level tended to spread the 
energy from the jet over the entire model grid and added erroneous structure in areas 
where the temperature gradients did not support it. The data weights (?/) were incre- 
mentally varied for each individual case. A weight was then chosen for each case to 
optimize a solution in the energy range 200-300 km 2 d -2 , in order to achieve solutions 
that were the most qualitatively correct with respect to the imagery and the CTZ88 Grid 
II ADCP velocity data. The values for y varied from .03-.49 with the weights in 15 of 
the 23 cases falling in the range .03-. 06. 

Tables 4-8 show the characteristics of the solutions resulting from the addition of 
various combinations of the lines of synthetic velocity data as additional constraints on 
the inversion. 



Table 4. CHARACTERISTICS OF I98I2- 19823 WITH VELOCITY CON- 
STRAINTS 



Data Lines 
Added 


Energy 

(km 2 d 2 ) 


Divergence 

(d"‘) 


T„ (%) 


Krr (%) 


none 


180.0 


0.316 


5.94 


92.4 


E 


255.0 


0.389 


6.10 


92.1 


DE 


253.0 


0.376 


6.05 


88.6 


DEF 


271.0 


0.362 


6.05 


87.9 


EFG 


253.0 


0.345 


6.00 


89.0 


D-G 


271.0 


0.330 


6.01 


86.2 



There was no significant change in the percent temperature error, which in some 
cases actually got slightly larger by 1-2%. This was not entirely unexpected since the 
velocity constraints could act to produce a solution different from the most 
thermodynamically correct one. The velocity error was large for all solutions (80-100%) 
but was improved by 10-15% in some cases by the addition of the known velocity con- 
straints (image pairs 19812-19823, 19912-19923, 19812- 19912 and 19823-19923). 
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Table 5. CHARACTERISTICS OF 19823-19912 WITH VELOCITY CON- 



STRAINTS 



Data Lines 
Added 


Energy 

(km 2 d 2 ) 


Divergence 

(O 


T.„ (%) 


K,r (%) 


none 


259.0 


0.421 


34.20 


107.0 


F 


311.0 


0.494 


34.80 


110.0 


EF 


318.0 


0.433 


34.20 


106.0 



Table 6. CHARACTERISTICS OF 19912-19923 WITH VELOCITY CON- 



STRAINTS 



Data Lines 
Added 


Energy 

(km 2 d 2 ) 


Divergence 

(d-) 


T.„ (%) 


Krr (%) 


none 


293.0 


0.579 


7.76 


101.0 


F 


383.0 


0.654 


7.91 


103.0 


FG 


380.0 


0.612 


7.82 


102.0 



Table 7. CHARACTERISTICS OF 19812-19912 WITH VELOCITY CON- 



STRAINTS 



Data Lines 
Added 


Energy 

(km 2 d 2 ) 


Divergence 

(d- 1 ) 


T„ r (%) 


v m (%) 


none 


107.0 


0.314 


11.2 


97.9 


E 


244.0 


0.454 


13.8 


102.0 


EF 


244.0 


0.420 


13.1 


104.0 


DEF 


236.0 


0.407 


12.7 


98.3 


EFG 


243.0 


0.380 


12.4 


96.9 


C-F 


234.0 


0.417 


12.6 


96.9 


D-G 


256.0 


0.368 


12.4 


94.3 


B-G 


229.0 


0.347 


12.6 


90.4 
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Table 8. CHARACTERISTICS OF 19823-19923 WITH VELOCITY CON- 



STRAINTS 



Data Lines 
Added 


Energy 

(km 2 d 2 ) 


Divergence 

(d- 1 ) 


T,„ (%) 


Krr (%) 


none 


96.6 


0.276 


7.39 


96.0 


E 


232.0 


0.430 


8.87 


98.3 


EF 


224.0 


0.413 


8.36 


98.7 


DEF 


226.0 


0.389 


8.42 


90.3 


EFG 


228.0 


0.354 


8.03 


90.5 


C-F 


231.0 


0.416 


8.47 


89.9 


D-G 


247.0 


0.332 


8.35 


86.2 


B-G 


232.0 


0.330 


8.83 


83.5 



Adding the velocity data greatly improved the ability of the model to reproduce 
mesoscale features that were apparent in the in situ data. This was demonstrated by the 
data addition sequence for image pair 19812-19823 in which velocity data was cumula- 
tively added to the inverse model (Figures 12-16). Data line E was generated from the 
ADCP data most coincident to this image pair with data lines D and F the next most 
coincident. The addition of just one line of synthetic velocity data significantly improved 
the definition of the jet by eliminating the "hole" in the middle of the inversion-alone 
solution (Figure 12). Steady improvement resulted from the addition of more data 
(Figures 13-16) but the biggest qualitative improvement came through the addition of 
just line E. Figures 12 and 16 look almost the same. The model jet was wider 
(55 vs 74 km) than the jet observed in the ADCP data (Figure 2) and located a bit 
further to the north. The model also failed to produce any indication of the return flow 
to the south of the jet and had some erroneous nearshore velocities to the south of Point 
Arena. The missing return flow was likely due to a lack of strong gradients in the 
warmer water of the return flow region, and the nearshore problems could be attribut- 
able to coastal fog. The best percent misfit to the actual velocity data (83%) was 
achieved by using all of the synthetic velocity data with image pair 19823-19923 
(Figure 17). Comparing this with the best results from image pair 19812-19923 
(Figure 16) showed the two solutions to be qualitatively similar, but differing in detail. 
Figure 16 produced a better representation of the anticyclonic eddy in the northwest, 
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while Figure 17 had a more coherent convergent flow leading into the jet from the 
north. Although Figure 17 could appear to have too much energy, the fact that velocity 
comparisons were done only at points where ADCP data existed would allow that sol- 
ution to be quantitatively the best but perhaps be qualitatively lacking. The solution 
from this image pair resulting from the addition of data lines D, E and F was considered 
to be the most qualitatively correct (Figure 18). The jet in the vicinity of the cold fila- 
ment was well defined but more narrow than in some of the other solutions. Addi- 
tionally, spurious structure in the low gradient areas of the images was kept to a 
minimum. 

Relatively poor results are depicted in Figures 19 and 20. The inversion-alone sol- 
ution in this case was the worst of the five cases, having the largest temperature and 
velocity errors. Possible reasons for the large errors were discussed in the previous sec- 
tion. Predictably, addition of the synthetic velocity data to the problem did little to im- 
prove it. 

A drawback to this technique was that it only constrained the velocity component 
normal to the synthetic GEOSAT track. Descending GEOSAT orbits cannot be used 
near the western coast of the United States owing to data problems as the satellite 
transitioned from over land to over water. For this and other similar coastal areas, de- 
scending orbit data that would help define the radar altimeter derived current velocities 
would not be available and only one velocity component could be determined. Fortu- 
itously, this synthetic data was nearly normal to the isotherms and the jet axis thereby 
providing the most needed (along-isotherm) information. 

Another slighly unrealistic aspect of the synthetic data was its physical and temporal 
density. The spacing of the tracks in this project was 25 km as compared to the 
= 130 km spacing of GEOSAT ground tracks. There w’ould also have been some 
temporal spacing to the GEOSAT ground tracks. This means that only one line of ac- 
tual data would be really useful in combination with the inverse model and its utility 
would be directly related to the temporal separation between the GEOSAT data and the 
imagery used with the inverse model. Use of more than one track of GEOSAT data 
would requre additional temporal weighting on the second and subsequent tracks that 
would likely cause them to have negligible effect on the inverse model. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



An inverse model using the heat equation with AVHRR imagery and constraints 
on the divergence and energy of the solution and with additional constraints in the form 
of synthetic GEOSAT-derived geostrophic current velocities showed promise as a means 
of improving the determination of near-surface current velocities entirely by remote 
means. In most of the cases run, both qualitative (visual comparison of current vectors) 
and quantitative (percent misfit to known velocities) improvement in the model's veloc- 
ity solutions was obtained when synthetic velocity data was added. The model per- 
formed best with an image separation of 12 hours which minimized the translation of 
features over the interval between images. Excessive feature translation between the 
images smeared the temperature gradients resulting in incorrect velocities in the model's 
solutions. The physical effects of the addition of the synthetic velocity data did not 
spread much beyond the actual points at which it was added. For this reason, the 
qualitative improvement in some cases appeared to be better than was indicated by the 
quantitative measurements. There was also a tendency for the synthetic velocity data 
to misrepresent the jet since the ADCP data was collected over a six day period during 
which the jet location was changing. 

The narrowness and orientation of the jet with respect to the model grid appeared 
to cause "holes" in the model's solutions. The model's apparent problem in these cases 
with the small width of the cold filament could be resolved by changing the model's tile 
size. With an image interval of 12 hours, 8x8 pixel tiles could be used without violating 
the CFL criteria discussed in Chapter I. Halving the tile size would double the size of 
the coefficient matrix with a concommitent increase in the computer processing time for 
the model. 

The model could be improved for this data set by experimenting with the weights 
on the divergence and energy constraints (a and /?) in order to find a set that improved 
the inversion-alone solutions. This is a process that would lend itself well to automation. 
Iterative computer techniques could be used to select weights for the constraints based 
on misfit levels chosen as a result of error analysis performed on the temperature data 
and the radar altimeter derived velocity data. It remains to be seen whether divergence 
and energy weights must be chosen for each image pair or u’hether one set of weights 
can apply to a set of images closely spaced in time or even closely related through similar 
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gradients. Using the inverse model with synthetic data generated by advecting a real 
image with a more sophisticated numerical model than that used by Kelly (1989) may 
help refine the selection criteria for the weights on the inversion constraints. 

Another possible improvement would be to weight the velocity data constraints 
based on the sine of the angle between the satellite track and the temperature gradients. 
This would allow only the model to define the cross-isotherm velocities when the radar 
altimeter track was parallel to the temperature gradients. Conversely, when the track 
was normal to the temperature gradients, only the radar altimeter derived velocities 
would drive the along-isotherm velocities which are in the model's null space. 

In order to incorporate real radar altimeter-derived geostrophic velocities into the 
model, a temporal weighting scheme would have to be used on the velocity data, such 
that radar altimeter data most coincident in time with the AVHRR data are weighted 
more than the less synoptic data. 

The value of an operational system based on these methods increases in direct pro- 
portion to the amount and accuracy of data available. Future refinements to the geoid 
will increase the accuracy of the geostrophic velocities obtainable from radar altimeter 
data. Follow-on radar altimeters to GEOSAT may have better temporal and physical 
track separation which would result in more altimeter data coincident to any set of se- 
quential infrared imagery. More than one altimeter in orbit would increase the accuracy 
of such a system even more by providing multiple coincident altimeter tracks per 
AVHRR image. Closely coincident tracks from ascending and descending orbits would 
be nearly normal to one another, allowing computation of velocity vectors rather than 
just the one normal component used in this project. 
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APPENDIX A. FIGURES 




Center of 
Voss 



Fig. I. Satellite altimetry relationships: Satellite height h subtracted from the 

satellite orbit height r yields the variation of sea surface height relative to 
the center of the earth. The reference ellipsoid is the best smooth approx- 
imation to the gcoid. [After Stewart, (1985)| 




Fig. 2. ADCF velocity field from CTZ88 Grid II: [Alter lluyer ct al., (1990)] 
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Longitude 



Fig. 3. Dynamic height field from CTZ88 Grid II: [After Iluycr ct al., (1990)] 



24 




of the backscattcrcd signal of the four beams measures the relative velocity 
between the ship and the ocean. Range-gating the signal allows determi- 
nation of currents at different depths. [Al ter Kosro, (19S5)] 
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Fig. 5. Synthetic GEOSAT-derived velocity field: Line B is to the northwest 

proceeding to G in the southwest 
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Fig. 6. CTZ88 Grid II ADCP Field on Image 19812: Line A is to the northwest 

proceeding to G in the southwest 
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Fig. 8. 19823-19912 Inversion-alone solution on Image 19912 
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Fig. 9. 19912-19923 Inversion-alone solution on Image 19923 
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Fig. 11. 19823-19923 Inversion-alone solution on Image 19912 
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Fig. 13. 19812-19823 Solution with data lines DE added 
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Fig. 14. 19812-19823 Solution with data lines DEF added 
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Fig. 15. 19812-19823 Solution with data lines EFG added 
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Fig, 17. 19823-19923 Solution with all velocity data added 
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Fig. 18. 19812-19823 plus DEF data lines on Image 19823 
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